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Abstract 



Analytical expressions correlating the volumetric flow rate to the inlet and 
outlet pressures are derived for the time-independent flow of Newtonian fluids in 
cylindrically-shaped elastic tubes using a one-dimensional Navier-Stokes flow model 
with two pressure-area constitutive relations. These expressions for elastic tubes 
are the equivalent of Poiseuille and Poiseuille-type expressions for rigid tubes which 
were previously derived for the flow of Newtonian and non-Newtonian fluids under 
various flow conditions. Formulae and procedures for identifying the pressure fleld 
and tube geometric proflle are also presented. The results are validated by a flnite 
element method implementation. Sensible trends in the analytical and numerical 
results are observed and documented. 

Keywords: fluid mechanics; Navier-Stokes; one-dimensional flow; Newtonian 
fluids; cylindrical elastic tubes; flnite element; time-independent; blood flow. 
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1 Introduction 

Considerable amount of work has been done in the past on the flow in rigid tubes 
with different types of geometry for both Newtonian and non-Newtonian fluids 
using various derivation methods (see for example [1-8]). However, relatively little 
work has been done on the flow in elastic tubes especially on developing closed- 
form analytical relations. These relations are useful in many scientific, industrial 
and medical applications; an obvious example is the flow of blood in large vessels. 
Most of the reported work in the literature on the flow in elastic tubes is based on 
the use of numerical methods such as finite element (see for instance [9, 10]) mainly 
due to the fact that since the flow in networks of elastic tubes was the main focus 
of these studies numerical methods were more appropriate to use. 

In the current paper, explicit analytical relations linking the volumetric flow 
rate to the pressure at the inlet and outlet are derived from a one- dimensional form 
of the Navier-Stokes equations for cylindrically-shaped elastic tubes with constant 
cross sectional area using two pressure-area constitutive models. The flow rate 
formulae are validated by a finite element implementation based on a Galerkin 
method with Lagrange polynomial interpolation and Gauss quadrature integration 
schemes. Formulae implicitly defining the tube profile and pressure field at each 
point along the tube axis are also provided, demonstrated and validated. The 
results presented in this paper are especially useful in biological studies such as 
modeling blood flow in arteries and veins. 

2 One- Dimensional Navier-Stokes Flow Model 

The widely used one- dimensional Navier-Stokes model describing the flow of New- 
tonian fluids, which is mainly formulated to model the flow in elastic tubes, is given 
by the following mass and momentum conservation principles 
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dt dz \ A J p dz A 

where A stands for the tube cross sectional area, t for time, Q for the volumetric 
flow rate, z for the space coordinate along the tube axis, L for the length of tube, 
a (= -^ " 2 with u and u being the fluid local and mean axial speed at the tube 
cross section respectively) for the axial momentum flux correction factor, p for the 
fluid mass density, p for the 2;-dependent pressure, and k for the viscosity friction 
coefficient which is usually given by k = 27raz//(a — 1) with u being the fluid 
kinematic viscosity defined as the ratio of the dynamic viscosity p to the mass 
density [9-13]. 

In using this model we assume a laminar, axi-symmetric, Newtonian, incom- 
pressible, fully-developed flow with insignificant gravitational forces and no-slip- 
at-wall boundary conditions [13, 14]. For our current purpose, we also assume a 
steady time-independent fiow and hence we drop the time terms in the Navier- 
Stokes equations. In the following section, we use this one-dimensional Navier- 
Stokes formulation to derive Q-p relations for cylindrical elastic tubes using two 
pressure-area constitutive relations. 

3 Deriving Q-p Relations 

For time independent fiow, the Navier-Stokes system given by Equations 1 and 2, 
becomes 
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dQ 

dz 
d faQ^\ Adp Q 
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The first of tliese equations states tliat Q as a function of z is constant. With 
regard to the second equation we have 

Adp _AdpdA _ d f A dp dA _ d f A dp 
p dz p dA dz dz J p dA dz dz J p dA 

Hence Equation 4 becomes 

d faQ^ f Adp ^\ Q 
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3.1 First p-A Model 

For this p-A model we assume a hnear pressure-area constitutive relation and hence 
the pressure is proportional to the change in cross sectional area relative to the 
reference area, that is 

p = 7(A-A„) (7) 

where p is the actual pressure as opposed to the reference pressure to which the 
reference area is defined, 7 is the proportionality coefficient which correlates to 
the tube stiffness, A is the tube cross sectional area at pressure p, and Ao is the 
reference area as identified by the reference pressure which, in this equation, is 
set to zero for convenience without affecting the generality of the results. From 
Equation 7, we get g^ = 7 and therefore 
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where the constant of integration is neglected because it will eventually vanish 
by the action of z partial derivative operator in Equation 6. Hence Equation 6 
becomes 
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that is 
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dz aQ 7^2 
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dA kA npQ 

On integrating both sides of this equation with respect to A, which is always 
greater than zero, we get 

. = H«l„^-_I^^» + C (13) 

K 6KpQ 

where C is the constant of integration which can be obtained from one of the two 
boundary conditions, e.g. the inlet boundary condition where A = Ain at z = 
with Ain being the tube inlet area, that is 



^-v'^-^'-^s^-^?" ("' 
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Hence 



Now, from this equation combined with the other boundary condition which 
defines the pressure at the outlet, that is A = A^u &X z = L where Ao„ is the tube 
outlet area and L is the tube length, we obtain 



This equation can be transformed to a quadratic polynomial in Q, i.e. 



In (^) Q' - ^Q + ^ (^i - ^L) = (17) 



with the following two roots 

_L±^L^-^l In (A,„M„0 3^5 (A3„ - AU 

^ ~ 2fln(^„/An) ^^^^ 

For Ain > Aou, which can always be satisfied by proper labeling, the two roots 

are necessarily real. For a physically viable fiow consistent in direction with the 

pressure gradient the negative sign should be taken to obtain a positive fiow and 

hence 

^ L _ ^^2 _ 4^ In (A,„M„0 jL. (^Al - A3J 
^ ^ 2fln(^„/An) ^^^^ 

This is due to the fact that for A^ > Aou the denominator is negative and hence 
to obtain a positive fiow rate the numerator should be negative as well, which is 
the case only if the negative sign is taken because the square root is always greater 
than L. This relation for elastic tubes is the equivalent of the Poiseuille equation 
for rigid tubes. However, for elastic tubes the fiow rate is dependent not on the 
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pressure difference but on the actual pressure at the inlet and outlet. 

3.2 Second p-A Model 

For the second pressure-area constitutive relation, the pressure is proportional to 
the radius change with a proportionality stiffness factor scaled by the reference 
area, that is 

P=£(v^-V^) (20) 

where p is the pressure, /3 is the tube stiffness factor, A^ is the reference area at 
the reference pressure and A is the area at pressure p. The tube stiffness factor for 
the second p-A model is normally defined by the following relation 



TxhoE 
1^? 



/3 = ^^ (21) 



where ho is the tube wall thickness at reference pressure, and E and q are respec- 
tively the Young's elastic modulus and Poisson's ratio of the tube wall. 

From the pressure-area constitutive relation of Equation 20 we obtain ^ = 



„ , n> and therefore 



Adp ^ ^ [A 13 ,, P ,o/2 , , 
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where the constant of integration is ignored as in the past. Hence Equation 6 
becomes 



dz 
that is 
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d faQ^ (3 ./2\ dA Q 



I.e. 



dA\ A 3pAo J dz A 



o:Q^ + ^^vA ^ + .^ . (25) 



A2 2pAo J dz A 
Following similar steps to those outlined in the first model, we obtain 

aQHniA/A,^)-^(A^/'-A%') 
z = ^ ^ (26) 

From the last equation associated with the second boundary condition at the 
outlet, i.e. A = Aou at z = L, we obtain the following expression for the volumetric 
flow rate 



-kL ± Jk^L^ - 4a In {A,JA,^) ^^ (^AIL' - A^i') 
^ ^ 2aln{A,JAou) ^^^^ 

Both these solutions are necessarily real for Am > Aou which can always be sat- 
isfied for normal flow conditions by proper labeling. For a flow which is physically- 
consistent in direction with the pressure gradient, the root with the plus sign should 
be selected, i.e. 



-kL + Jk^L^ - 4«ln (A^JAou) ^^ (aIL' - A^i' 

Q y '_ /OS 

^ 2aln(A,„/A,) ^ 

This, in essence, is a relation between flow rate and pressure drop, similar to 
the Poiseuille law for rigid tubes, although for elastic tubes the flow rate, as given 
by Equation 28, does not depend on the pressure difference, as for rigid tubes, 
but on the actual inlet and outlet pressure as defined by the inlet and outlet area 
respectively. 
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4 Finite Element Formulation 

The flow formulae derived in the previous section can be validated by the finite 
element method using the weak formulation. This formulation is outlined for the 
first and second p-A models in the following two subsections. More details about 
the finite element technicalities and the solution scheme using Newton-Raphson 
iteration are given in [13]. 

4.1 First p-A Model 

The Navier-Stokes system, given by Equations 1 and 2, can be cast in matrix 
form which is more appropriate for numerical manipulation and implementation as 
follow 



dV OF 
dt dz 



B = 



(29) 



where 



U 



A 

Q 



Q 






and 



B 



K- 



(30) 



On multiplying Equation 29 by weight functions and integrating over the solu- 
tion domain, z, the following system is obtained 



(9U 



dF 



■ u)dz + / — — ■ u}dz + / B • u}dz = 
n dt Jn dz J^ 



(31) 



where Q is the solution domain, and c*; is a vector of arbitrary test functions. On 
integrating the second term of Equation 31 by parts, the following weak form of 
the preceding ID flow system is obtained 
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/ — - • ujdz - / F ■ -^dz + / B ■ ujdz + [F ■ uj]on = 
Jn dt J^ dz J^ 



(32) 



where dQ is the boundary of the solution domain. This weak formulation, cou- 
pled with suitable boundary conditions, can be used as a basis for finite element 
implementation in conjunction with an iterative scheme such as Newton-Raphson 
method. Following a solution scheme detailed in [13] and based on the method 
of characteristics [9, 10, 15-18], the eigenvalues Ai,2 and left eigenvectors Li^2 of 
the H matrix, which are required for obtaining the compatibility conditions on the 
boundaries, are found as follow 



det (H - AI) = det 



/ 



VL 



A2 



+ 



■yA 2aQ 



X 



(33) 



where H is the matrix of partial derivatives of F with respect to U, that is 



H 



OF 
dV 



■ A2 



+ 



^A 2aQ 



(34) 



On solving Equation 33 the eigenvalues are obtained 



Ai, 



aQ Q\ ^ jA 



(35) 



which are necessarily real for a > 1 as it is always the case, and hence the left 
eigenvectors are obtained 



il, 



-a2± 



A2 



(a2-a) + 2^ 1 



■yA 



(36) 



The compatibility conditions for the time-independent flow arising from pro- 
jecting the differential equations in the direction of the outgoing characteristic 
variables at the inlet and outlet are then obtained from 
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Xi,2(H^ + B) =0 



that is 
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Q 

'A 
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4.2 Second p-A Model 

Following a similar procedure to that outlined in the previous subsection for the 
first p-A model, the finite element formulation leads to the following matrix struc- 
ture, eigenvalues, left eigenvectors and time-independent compatibility conditions 
respectively 



U 



A 

Q 



Q 



A "•" 3pAo 



and 



B 



K, 



Q 



(40) 
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a— ± 
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and 



Q , Q\ 2 .^PVA\dQ Q\pVA\dAfdQ\Q 

-a— ± \ -— (a'' — a) -^ :— -^ — h — «-— - H — — — h 2a— — h k — = 

A \l A^^ ' 2pAo dz y A^ 2pAoj dz \ dz ) A 

(43) 



5 Numerical Validation 

To validate the derived flow formulae, the finite element formulation as outlined 
in the previous section was implemented for the two p-A models in a computer 
code using a Galerkin method with a Lagrange polynomial interpolation associated 
with a Gauss quadrature integration scheme. The comparison between the analytic 
and finite element solutions is outlined for some typical cases in the following two 
subsections. 

5.1 First p-A Model 

Extensive tests have been carried out to verify Equation 19; a sample of which is 
given in Table 1. Certain sensible trends can be observed in these results. For 
example, the diagonally-oriented entries from top-left to bottom-right direction in 
the table are of similar magnitude which is sensible since in this quasi-linear flow 
regime obtained at relatively low pressures the flow is Poiseuille-like and hence it 
is almost proportional to the pressure difference (i.e. Pin — Pou)- This Poiseuille- 
like behavior disappears at high-pressure flow regimes as the flow rate becomes 
increasingly dependent on the actual pressures at the inlet and outlet rather than 
on the pressure difference. Another sensible trend is that the flow rate in these 
diagonally-oriented entries is increasing in the top-left to bottom-right direction 
due to the fact that although the pressure difference for these entries is the same. 
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the lower entries have larger area at the inlet and outlet, due to the higher pressure 
at the tube entrance and exit, than the upper ones. This trend is more obvious at 
higher pressure regimes. 

We also used Equation 15, which implicitly correlates A to z, to obtain the 
pressure field inside the tube and the tube profile by numerically solving for A 
for a given z. A sample of these results, with their finite element counterparts, is 
presented in Figures 1 and 2. These figures confirm the sensibility of the obtained 
analytical and numerical results. 



Table 1: Sample results of the volumetric fiow rate in m^/s related to the elastic 
tube investigation for the first p-A model. The rows stand for the inlet pressure. 
Pin, and the columns for the outlet pressure, Pou, in Pa- The parameters with 
which these results are obtained are: p = 1060 kg/m^, /i = 0.0035 Pa.s, a = 1.333, 
L = 1.0 m, r = 0.1 m, and 7 = 5 x 10^ Pa/m^. In each Pj„ row the top and 
bottom entries are respectively the analytic solution, given by Equation 19, and 
the finite element solution which is obtained with a quadratic Lagrange polynomial 
interpolation. 



Pi,t 





100 


200 


300 


Po 
400 


500 


600 


700 


800 


900 


100 


0.286046 
0.286046 


200 


0.307977 
0.307977 


0.286332 
0.286332 


















300 


0.315789 
0.315789 


0.308278 
0.308278 


0.286619 
0.286619 
















400 


0.319850 
0.319850 


0.316096 
0.316096 


0.308579 
0.308579 


0.286905 
0.286905 














500 


0.322373 
0.322373 


0.320158 
0.320159 


0.316402 
0.316402 


0.308881 
0.308881 


0.287192 
0.287192 












600 


0.324118 
0.324119 


0.322684 
0.322684 


0.320467 
0.320467 


0.316708 
0.316708 


0.309182 
0.309182 


0.287479 
0.287479 










roo 


0.325415 
0.325415 


0.324430 
0.324430 


0.322994 
0.322994 


0.320776 
0.320777 


0.317015 
0.317015 


0.309484 
0.309484 


0.287766 
0.287766 








800 


0.326430 
0.326430 


0.325727 
0.325728 


0.324741 
0.324742 


0.323305 
0.323305 


0.321086 
0.321086 


0.317322 
0.317322 


0.309785 
0.309785 


0.288053 
0.288053 






900 


0.327256 
0.327257 


0.326743 
0.326743 


0.326040 
0.326040 


0.325053 
0.325053 


0.323615 
0.323616 


0.321396 
0.321395 


0.317628 
0.317628 


0.310087 
0.310087 


0.288340 
0.288340 




1000 


0.327950 
0.327951 


0.327569 
0.327570 


0.327056 
0.327056 


0.326352 
0.326353 


0.325365 
0.325365 


0.323926 
0.323926 


0.321704 
0.321704 


0.317935 
0.317935 


0.310389 
0.310389 


0.288627 
0.288627 
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Figure 1: Pressure versus tube axial coordinate for three sample cases related to 
the first p-A model as obtained analytically from Equation 15 and numerically by a 
finite element method with a quadratic polynomial interpolation scheme as outlined 
in the previous section. The labels '1', '2' and '3' in these plots refer respectively 
to the cases where P^ = 1000 Pa and Pou = Pa, Pj„ = 500 Pa and Pou = Pa, 
and Pin = 700 Pa and Pqu = 400 Pa. The tube, fluid and flow parameters with 
which these results are obtained are: p = 1060 kg/m^, /x = 0.0035 Pa.s, a = 1.333, 
L = 1.0 m, r = 0.1 m, and 7 = 5 x 10^ Pa/m^. 
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Figure 2: Radius versus tube axial coordinate for the three sample cases of Figure 
1. 
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5.2 Second p-A Model 

Extensive tests have been carried out to verify Equation 28; a sample of which is 
given in Table 2. Also, we used Equation 26 to obtain the pressure field inside the 
tube and the tube geometric profile, as outlined for the first p-A model. A sample 
of these results, with their finite element equivalents, is presented in Figures 3 and 
4. Similar sensible trends to those observed in the first p-A model are detected. 



Table 2: Sample results of the volumetric flow rate in m^/s related to the elastic 
tube investigation for the second p-A model. The rows stand for the inlet pressure. 
Pin, and the columns for the outlet pressure, Pou, in Pa. The parameters with 
which these results are obtained are: p = 1060 kg/m^, /i = 0.0035 Pa.s, a = 1.333, 
L = 1.0 m, r = 0.1 m, and /3 = 5 x 10'^ Pa.m. In each Pj„ row the top and 
bottom entries are respectively the analytic solution, given by Equation 28, and 
the finite element solution which is obtained with a quadratic Lagrange polynomial 
interpolation. 



Pin 





100 


200 


300 


Po 
400 


500 


600 


700 


800 


900 


100 


0.273135 
0.273135 


200 


0.292950 
0.292950 


0.273397 
0.273397 


















300 


0.299986 
0.299986 


0.293221 
0.293221 


0.273659 
0.273659 
















400 


0.303637 
0.303637 


0.300259 
0.300259 


0.293491 
0.293491 


0.273922 
0.273922 














500 


0.305904 
0.305904 


0.303912 
0.303912 


0.300532 
0.300533 


0.293762 
0.293762 


0.274184 
0.274184 












600 


0.307471 
0.307471 


0.306180 
0.306180 


0.304187 
0.304187 


0.300806 
0.300806 


0.294033 
0.294033 


0.274447 
0.274447 










700 


0.308634 
0.308634 


0.307747 
0.307747 


0.306455 
0.306456 


0.304462 
0.304462 


0.301080 
0.301080 


0.294304 
0.294304 


0.274710 
0.274710 








800 


0.309543 
0.309543 


0.308910 
0.308910 


0.308023 
0.308023 


0.306731 
0.306731 


0.304737 
0.304737 


0.301353 
0.301353 


0.294575 
0.294575 


0.274973 
0.274973 






900 


0.310283 
0.310283 


0.309820 
0.309820 


0.309187 
0.309187 


0.308299 
0.308300 


0.307007 
0.307007 


0.305012 
0.305012 


0.301627 
0.301627 


0.294847 
0.294847 


0.275236 
0.275237 




1000 


0.310903 
0.310904 


0.310560 
0.310561 


0.310097 
0.310097 


0.309464 
0.309464 


0.308576 
0.308576 


0.307283 
0.307283 


0.305287 
0.305287 


0.301902 
0.301902 


0.295118 
0.295118 


0.275500 
0.275500 
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Figure 3: Pressure versus tube axial coordinate for three sample cases related to the 
second p-A model as obtained analytically from Equation 28 and numerically by a 
finite element method with a quadratic polynomial interpolation scheme as outlined 
in the previous section. The labels '1', '2' and '3' in these plots refer respectively 
to the cases where Pm = 900 Pa and Pou = Pa, Pm = 400 Pa and Pou = Pa, 
and Pin = 600 Pa and Pou = 200 Pa. The tube, fluid and flow parameters with 
which these results are obtained are: p = 1060 kg/m^, /x = 0.0035 Pa.s, a = 1.333, 
L = 1.0 m, r = 0.1 m, and /3 = 5 x 10^ Pa.m. 
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Figure 4: Radius versus tube axial coordinate for the three sample cases of Figure 
3. 
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6 Conclusions 

In this study, two analytical expressions, correlating volumetric flow rate to pres- 
sure at inlet and outlet, are derived for the Newtonian flow in cylindrical elastic 
tubes from a one-dimensional form of the Navier-Stokes equations. The expressions 
are validated numerically by a finite element method based on a Galerkin scheme 
with Lagrange interpolation and Gauss quadrature integration. Sample results, 
which are quantitatively and qualitatively sensible, are presented for demonstra- 
tion. Two constitutive relations, depicting the nature of the relation between area 
and pressure in elastic tubes, are used in all these derivations and finite element 
implementation. The foundations of the finite element weak form for the two p-A 
models are outlined for completion. Preliminary rational trends in these results 
are observed and documented. Analytical implicit relations for obtaining the pres- 
sure field inside the tube, as well as the tube geometric profile, are also presented, 
demonstrated and numerically validated. The outcome of this investigation, nu- 
merical as well as analytical, is of relevance to several areas of science, technology 
and medicine. 
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Nomenclature 

a correction factor for axial momentum flux 

P stiffness factor in the second p-A model 

7 stiffness factor in the first p-A model 

K viscosity friction coefficient 

A 1,2 eigenvalues of H matrix 

H fluid dynamic viscosity 

u fluid kinematic viscosity 

p fluid mass density 

q Poisson's ratio of tube wall 

u} vector of test functions in flnite element formulation 

Q solution domain 

dQ boundary of solution domain 

A tube cross sectional area at pressure p 

Ain tube cross sectional area at inlet 

Ao tube reference cross sectional area at reference pressure 

Aou tube cross sectional area at outlet 

B matrix of force terms in the ID Navier-Stokes equations 

E Young's modulus of tube wall 

F flux matrix in the ID Navier-Stokes equations 

H matrix of partial derivatives of F with respect to U 

ho tube wall thickness at reference pressure 

L length of tube 
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-Li 2 left eigenvectors of H matrix 

p pressure at given coordinate z 

Pin pressure at tube inlet 

Pou pressure at tube outlet 

Q volumetric flow rate 

r radius 

t time 

u local axial speed of fluid at cross section 

u mean axial speed of fluid at cross section 

U vector of Navier-Stokes dependent variables 

z tube axial coordinate 
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